#### LOAD LIBRARIES
source("/home/jose/Documents/funciones2.R")
packages <- c("dplyr", "chron", "ggplot2", "ggdist", "emmeans", "ggsci", "DHARMa", "glmmTMB", "performance")
lapply(packages, library, character.only = TRUE)

######################################################################
##### 															                             #####
##### 															                             #####
##### **************** CTMAX AND HIGH TEMP EXPS **************** #####
##### 														                               #####
##### 														                               #####
######################################################################

##### LOAD DATA ####
ctmax_data <- read.csv("/home/jose/Documents/papers/mios/7-de-oro/para-paper/ctmax_data.csv", sep = ",", dec = ".")

#### FIT LINEAR MODEL ####
fit <- lm(t_coma ~ assay_temp * sex, data = ctmax_data)

#### TEST FOR SEX DIFFERENCES
anova(fit)

#### SINCE NO DIFFERENCES ON SEX, WE MERGE DATA
#### FIT LINEAR MODEL
fit <- lm(t_coma ~ assay_temp, data = ctmax_data)

#### TEST FOR SEX DIFFERENCES
anova(fit)

#### CALCULATE Z COEFFICIENT
-1 / coefficients(fit)[2] # Z 7.445772

#### CALCULATE CTMAX
-coefficients(fit)[1] / coefficients(fit)[2] # CTmax 70.3145

######################################################################
###################### PROPORTION OF RECOVERY ########################
######################################################################

high_temp <- read.csv("/home/jose/Documents/papers/mios/7-de-oro/para-paper/high_temp_data.csv", sep = ",", dec = ".")

#### CREATE MODEL FOR RECOVERY ####
m.recov.high <- glm(rec30 ~ sex * assay_temp, family = binomial, data = high_temp)

##### CHECK ASSUMPTIONS #####
sim <- simulateResiduals(m.recov.high, n = 1000)
plotQQunif(sim)
plotResiduals(sim)

testDispersion(sim)
testZeroInflation(sim)

#### OBTAIN ANOVA STATS ####
anova(m.recov.high, test = "Chisq")

#### CALCULATE TRENDS AND COMPARE SEXES ####
(comp.pend <- emtrends(m.recov.high, "sex", var = "assay_temp"))
pairs(comp.pend)


######################################################################
###################### PROPORTION OF SURVIVAL ########################
######################################################################

#### CREATE MODEL FOR SURVIVAL ####
m.sup.high <- glm(sup24 ~ sex * assay_temp, family = binomial, data = high_temp)

##### CHECK ASSUMPTIONS #####
sim <- simulateResiduals(m.sup.high, n = 1000)
plotQQunif(sim)
plotResiduals(sim)

testDispersion(sim)
testZeroInflation(sim)

#### OBTAIN ANOVA STATS ####
anova(m.sup.high, test = "Chisq")

#### CALCULATE TRENDS AND COMPARE SEXES ####
(comp.pend <- emtrends(m.sup.high, "sex", var = "assay_temp"))
pairs(comp.pend)


######################################################################
#####															                               #####
#####															                               #####
##### ********************* CHILL RECOVERY ********************* #####
#####														                                 #####
#####														                              	 #####
######################################################################


######################################################################
######################### TIME TO RECOVERY ###########################
######################################################################

chill <- read.csv("/home/jose/Documents/papers/mios/7-de-oro/para-paper/chill-data.csv", sep = ",", dec = ".", as.is = F)

#### CREATE MODELS ####
m.chill <- glm(recovery.time ~ time * sex, family = Gamma(link = "identity"), data = chill)
m.chill2 <- lm(log(recovery.time) ~ time * sex, data = chill)
m.chill3 <- lm(log10(recovery.time) ~ time * sex, data = chill)

#### COMPARE MODELS ####
AIC(m.chill, m.chill2, m.chill3) # we select m.chill3

##### CHECK ASSUMPTIONS #####
sim <- simulateResiduals(m.chill3, n = 1000)
plotQQunif(sim)
plotResiduals(sim)
testDispersion(sim)
testZeroInflation(sim)

#### OBTAIN ANOVA STATS ####
anova(m.chill3, test = "Chisq")

######################################################################
############################## SURVIVAL ##############################
######################################################################

m.sup.chill <- glm(sup.24 ~ sex * time, family = binomial, data = chill)

##### CHECK ASSUMPTIONS #####
sim <- simulateResiduals(m.sup.chill, n = 1000)
plotQQunif(sim)
plotResiduals(sim)

testDispersion(sim)
testZeroInflation(sim)

#### OBTAIN ANOVA STATS ####
anova(m.sup.chill, test = "Chisq")

(comp.pend <- emtrends(m.sup.chill, "sex", var = "time"))
#### SE AND CIs ARE WRONGLY CALCULATED DUE TO PERFECT DISCRIMINATION
#### TURNING TO BAYESIAN APROXIMATION FOLLOWING GELMAN 2008 

library(arm)
m.sup.chill.bayes <- bayesglm(sup.24 ~ sex * time, data = chill, family = "binomial")

##### CHECK ASSUMPTIONS #####
sim <- simulateResiduals(m.sup.chill.bayes, n = 1000)
plotQQunif(sim)
plotResiduals(sim)

testDispersion(sim)
testZeroInflation(sim)

#### OBTAIN ANOVA STATS ####
anova(m.sup.chill.bayes, test = "Chisq")

#### CALCULATE TRENDS AND COMPARE SEXES ####
(comp.pend.bayes <- emtrends(m.sup.chill.bayes, "sex", var = "time"))
pairs(comp.pend)
